% This is the dynare code for Technical default model
% By Jun Yu 2021.5.30

% Endogenous variables
var 
% 18 Endogenous variables
y   % output
ii   % Investment 
rl  % loan rate
rd  % risk free rate on deposit
rk  % return on capital
q   % price of purchasing capital
m   % Household's SDF
n   % Net worth
b   % debt
act_def % actual default prob
spread
TDprob   % technical default probability
re   % equity return
premiumH
c   % Consumption
lamda   % Marginal value of net worth
omega0bar  % the default cutoff
h    % the debt-to-networth ratio
dk   % Growth rate of aggregate capital
a    % Exogenous TFP shock
dh   % household deposit
de   % entrepreneur deposit
sigma0   % sigma0 shock
sigma1   % sigma1 shock
u   % Household utility value
eu  % Household's expected utility value
%% 7 auxiliary variables
dn   % Growth rate of net worth
dc   % Growth rate of consumption
di   % Growth rate of investment
dy   % Growth rate of output
db   % Growth rate of debt
dq   % Growth rate of capital price
g    % Capital adjustment cost
gap
dTDprob
;
  
% in total 25 endo variables  

% Exogenous variables
%varexo ea e_0 e_1;
% varexo ea e_s;
varexo ea;

% Parameters 

parameters 
ppsi    % Relative risk aversion
beta    % discount rate
gamma   % Intertempral elasticity of substitution
a1      % capital adjustment cost function parameter
xi      % capital adjustment cost function parameter
a2      % capital adjustment cost function parameter
alppha  % capital share
delta   % capital depriciation rate
rho_a   % Persistence of TFP shock
sigma_a % STD of TFP shock
eta     % Entrepreneur servival rate
vphi    % bank's verfication or mornitoring cost
chi     % fraction of redrawed networth
mu_0    % mean of omega0 shock
sigma_0ss % STD of omega0 shock
mu_1    % mean of omega1 shock
sigma_1ss % STD of omega1 shock
rho_s1  % Persistence of omega1 shock
sigma_s1 % STD of sigma1 shock
sig0_per % fraction of sigma0 shock
sigma_s  % shock size
rho_s0   % Persistence of omega0 shock
sigma_s0 % STD of sigma0 shock
F_omega0 % the technical default probability

%% The steady state value as initial value
yss    % output
iiss   % Investment
rlss   % loan rate
rdss   % risk free rate on deposit
rkss   % return on capital
qss    % price of purchasing capital
mss    % Household's SDF
nss    % Net worth
bss    % debt
uss    % Household utility value
euss   % Household's expected utility value
css    % Consumption
lamdass   % Marginal value of net worth
dkss   % Growth rate of aggregate capital
ass    % Exogenous TFP shock
dhss   % household deposit
dess   % entrepreneur deposit
omega0barss  % the default cutoff
hss    % the debt-to-networth ratio
ress
PD2
;        

load para_v1;

set_param_value('ppsi',ppsi)
set_param_value('beta',beta)
set_param_value('gamma',gamma)
set_param_value('a1',a1)
set_param_value('xi',xi)
set_param_value('a2',a2)
set_param_value('alppha',alppha)
set_param_value('delta',delta)
set_param_value('ass',log(Ass))
set_param_value('rho_a',rho_a)
set_param_value('rho_s0',rho_s0)
set_param_value('rho_s1',rho_s1)
set_param_value('sigma_s1',sigma_s1)
set_param_value('sigma_s0',sigma_s0)
set_param_value('sigma_a',sigma_a)
set_param_value('eta',eta)
set_param_value('vphi',zeta)
set_param_value('chi',chi)
set_param_value('mu_0',mu_0)
set_param_value('sigma_0ss',log(sigma_0))
set_param_value('mu_1',mu_1)
set_param_value('sigma_1ss',log(sig1))
set_param_value('PD2',PD2)
set_param_value('sig0_per',sig0_per)
set_param_value('sigma_s',sigma_s)

%% steady state values
set_param_value('yss',log(Yss))
set_param_value('iiss',log(Iss))
set_param_value('rlss',log(RLss))
set_param_value('ress',log(REss))
set_param_value('rdss',log(RDss))
set_param_value('rkss',log(Rkss))
set_param_value('qss',log(Qss))
set_param_value('mss',log(Mss))
set_param_value('nss',log(Nss))
set_param_value('bss',log(Bss))
set_param_value('css',log(Css))
set_param_value('lamdass',log(lamda_ss))
set_param_value('dkss',log(DK))
set_param_value('ass',log(Ass))
set_param_value('dhss',log(Dhss))
set_param_value('dess',log(Dess))
set_param_value('omega0barss',omega0_barss)
set_param_value('hss',log(Hss))
set_param_value('uss',uss)
set_param_value('euss',euss)
set_param_value('F_omega0',F_omega0)

external_function(name = FUNC1, nargs = 6, first_deriv_provided, second_deriv_provided);
external_function(name = FUNC2, nargs = 6, first_deriv_provided, second_deriv_provided);
external_function(name = FUNC3, nargs = 6, first_deriv_provided, second_deriv_provided);
external_function(name = FUNC4, nargs = 6, first_deriv_provided, second_deriv_provided);
external_function(name = FUNC5, nargs = 6, first_deriv_provided, second_deriv_provided);

model;
%% Household Side
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% dc is actually the actual growth rate of unnormalized C
% dy, di, dj is of the same form
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define some auxilury variables
# Td_probm = normcdf(( omega0bar + 0.5*exp(2*sigma0)) / exp(sigma0));
# Td_probm1 = normcdf(( omega0bar(-1) + 0.5*exp(2*sigma0(-1))) / exp(sigma0(-1)));
# func6 = 1-normcdf(( log((exp(rl+h))/(exp(rk(+1))*(1+exp(h))))-omega0bar + 0.5*exp(2*sigma1)) / exp(sigma1));
# func7 = (1-normcdf(( log((exp(rl+h))/(exp(rk(+1))*(1+exp(h))))-omega0bar + 0.5*exp(2*sigma1)) / exp(sigma1)-exp(sigma1)));
# func8 = (normcdf(( log((exp(rl+h))/(exp(rk(+1))*(1+exp(h))))-omega0bar + 0.5*exp(2*sigma1)) / exp(sigma1)-exp(sigma1)));


% 1 HHs preference 
exp(u*(1-1/ppsi)) = (1-beta)*exp((1-1/ppsi)*c) + beta*exp(eu*(1-1/ppsi));

% 2 Auxiliary eu
exp((1-gamma)*eu) = exp((1-gamma)*(u(+1) + dk));

% 3 SDF m
m = log(beta) + (-1/ppsi)*(c - c(-1) + dk(-1)) + (1/ppsi - gamma)*(u - eu(-1) + dk(-1));
% m(+1) = log(beta) + (1-1/ppsi)*(cstar(+1) - cstar + dk) - (c(+1) - c + dk)  + (1/ppsi - gamma)*(u(+1) - eu + dk);

% 4 HH's Euler equation
1/exp(rd) = exp(m(+1));

% 5 Return on capital Rk
exp(rk) = (alppha* exp(a) + exp(q)*(1-delta))/exp(q(-1));

%% Capital good producer
% 6 capital adjustment cost
g = ln(a1/(1-1/xi)*exp((1-1/xi)*ii)+a2);

% 7 Capital price
q = 1/xi*ii- ln(a1);

% 8 Production Function
exp(y) = exp(a);

% 9 Law of motion of capital
exp(dk) = (1-delta) + exp(g);

% 10 TFP shock
a - ass = rho_a*(a(-1)-ass) - sigma_a*ea;

%% Entrepreneurs problem
% 11 Balance sheet condition
exp(n) + exp(b) = exp(q+dk) + exp(de);

% 12 Entrepreneur debt
exp(de)=Td_probm*(1+exp(h))*exp(n);

% 13 Debt
exp(b)=exp(n)*exp(h);

% 14 debt market clearing condition
exp(b) = exp(de) + exp(dh);

%% 15 Market Clearing condition
exp(y) = exp(c) + exp(ii)+ vphi*exp(n(-1)+rk)*(1+exp(h(-1)))*FUNC4(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)));

% 16 FOC of H 
exp(m(+1))*(exp(rk(+1))*FUNC1(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1)) -exp(rl)*FUNC3(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))) = exp(lamda)*(1-exp(rl+m(+1))*FUNC3(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))-(1-vphi)*exp(m(+1))*exp(rk(+1))*FUNC4(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))+vphi*exp(rl+m(+1))/(1+exp(h))*FUNC5(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))-exp(rd+m(+1))*Td_probm);

% 17 FOC of omega0bar
exp(m(+1))*(exp(rd)-exp(rk(+1))*(1+exp(h))*exp(omega0bar)*func7+exp(rl+h)*func6) = exp(lamda)*exp(m(+1))*(exp(rl+h)*func6+(1-vphi)*exp(rk(+1))*(1+exp(h))*exp(omega0bar)*func8-exp(rd+h));
%omega0bar = omega0barss;

% 18 FOC of RL
exp(m(+1))*(exp(h)*FUNC3(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))) = exp(lamda)*exp(m(+1))*(exp(h)*FUNC3(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))-vphi*exp(h)*FUNC5(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1)));

% 19 lender's break even condition
exp(h) = exp(m(+1))*(exp(rl+h)*FUNC3(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))+(1-vphi)*exp(rk(+1))*(1+exp(h))*FUNC4(omega0bar,exp(rl),exp(h),exp(rk(+1)),exp(sigma0),exp(sigma1))+exp(rd+h)*Td_probm);

% 20 Law of motion of net worth
exp(n)*exp(dk(-1)) = exp(n(-1))*eta*(exp(rk)*(1+exp(h(-1)))*FUNC1(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)))-exp(rl(-1)+h(-1))*FUNC3(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)))+exp(rd(-1))*Td_probm1)+chi*exp(n(-1))*(1-eta*(1-FUNC2(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)))));

% the actual default prob
act_def = FUNC2(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)));

% loan spread
spread = rl-rd;

% Auxiliary Growth rates 
% 21
dc = c - c(-1) + dk(-1);
% 22
dy = y - y(-1) + dk(-1);
% 23
di = ii - ii(-1) + dk(-1);
% 24
dn = n - n(-1) + dk(-1);
% 25
db = b - b(-1) + dk(-1);
% 30 
dq = q - q(-1);

% 31 return on equity
exp(re) = exp(rk)*FUNC1(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)))*(1+exp(h(-1)))-exp(rl(-1)+h(-1))*FUNC3(omega0bar(-1),exp(rl(-1)),exp(h(-1)),exp(rk),exp(sigma0(-1)),exp(sigma1(-1)))+exp(rd(-1))*Td_probm1;

% 32 Asset pricing equation
% exp(m(+1)+re(+1))=exp(lamda)*(1-Td_probm)+Td_probm;

% 33 External finance premium        
premiumH = rk(+1) - rd; % definiton of primium, Rk/R

%34 sigma0 Shock 
%log(sigma0/sigma_0) = rho_s0*log(sigma0(-1)/sigma_0)+sigma_s0*e_s;
% sigma0 - sigma_0ss = rho_s0*(sigma0(-1)-sigma_0ss) +sigma_s0*e_s;
%sigma0 - sigma_0ss = rho_s0*(sigma0(-1)-sigma_0ss) +sigma_s0*ea;
sigma0 = sigma_0ss;

%35 sigma1 Shock 
%log(sigma1/sigma_1)=rho_s1*log(sigma1(-1)/sigma_1)+sigma_s1*e_s;
%sigma1 - sigma_1ss = rho_s1*(sigma1(-1)-sigma_1ss) +sigma_s1*e_s;
%sigma1 - sigma_1ss = rho_s1*(sigma1(-1)-sigma_1ss) +sigma_s1*ea;
%sigma0 - sigma_0ss = rho_s0*(sigma0(-1)-sigma_0ss) +sigma_s0*ea;
sigma1 = sigma_1ss;

%36 the technical default probability
exp(TDprob) = normcdf(( omega0bar + 0.5*exp(2*sigma0)) / exp(sigma0));
%TDprob = normcdf(( omega0bar + 0.5*exp(2*sigma0)) / exp(sigma0));

dTDprob = TDprob - TDprob(-1);

gap = ( omega0bar + 0.5*exp(2*sigma0)) / exp(sigma0);

end;          
% in total 38 obs equations
steady_state_model;
%initval;
y =yss;  % output
ii = iiss;  % Investment
rl =rlss; % loan rate
rd =rdss; % risk free rate on deposit
re =ress; % equity return
rk =rkss; % return on capital
premiumH = rk-rd;
q  =qss; % price of purchasing capital
m  =mss; % Household's SDF
n =nss;  % Net worth
b =bss;  % debt
u =uss;  % Household utility value
eu =euss; % Household's expected utility value
c  =css; % Consumption
lamda =lamdass;  % Marginal value of net worth
dk =dkss;  % Growth rate of aggregate capital
a =ass; % Exogenous TFP shock
dh = dhss; % household deposit
de =dess;  % entrepreneur deposit
omega0bar = omega0barss; % the default cutoff
h =hss ; % the debt-to-networth ratio
g = iiss;
dc = dkss;
dy   = dkss;
di   = dkss;
dn   = dkss;
db = dkss;
dq = 0;
sigma0 = sigma_0ss;
sigma1 = sigma_1ss;
TDprob = log(F_omega0);
dTDprob = 0;
%TDprob = F_omega0;
act_def=PD2;
spread = rl-rd;
gap = ( omega0bar + 0.5*exp(2*sigma0)) / exp(sigma0);
end;

shocks;
var ea; stderr 1;
%var e_s; stderr 1; 
%var e_1; stderr 1; 
%corr e_s, ea = 0.8;
%corr e_1, ea = 0.99;
%corr e_0, e_1 = 1;
end;


steady;
check;
model_diagnostics;

write_latex_dynamic_model;
% set_dynare_seed(1000);
% options_.relative_irf
options_.pruning = 1;
stoch_simul(order=2,periods=2000,drop=800,irf=30, pruning,replic=1000,simul_replic=1000); 

%read out simulations
%simulated_series_raw=get_simul_replications(M_,options_);
%save('simulated_series_raw.mat','simulated_series_raw');

%filter series
%simulated_series_filtered=NaN(size(simulated_series_raw));
%for ii=1:options_.simul_replic
%    [trend, cycle]=sample_hp_filter(simulated_series_raw(:,:,ii)',1600);
%    simulated_series_filtered(:,:,ii)=cycle';
%end

%get variable positions
%y_pos=strmatch('y',M_.endo_names,'exact');
%c_pos=strmatch('c',M_.endo_names,'exact');
%i_pos=strmatch('invest',M_.endo_names,'exact');
%k_pos=strmatch('k',M_.endo_names,'exact');
%h_pos=strmatch('h',M_.endo_names,'exact');
%productivity_pos=strmatch('productivity',M_.endo_names,'exact');

%var_positions=[y_pos; c_pos; i_pos; k_pos; h_pos; productivity_pos];
%get variable names
%var_names=M_.endo_names_long(var_positions,:);

%Compute standard deviations
%std_mat=std(simulated_series_filtered(var_positions,:,:),0,2)*100;

%Compute correlations
%for ii=1:options_.simul_replic
%    corr_mat(1,ii)=corr(simulated_series_filtered(y_pos,:,ii)',simulated_series_filtered(y_pos,:,ii)');
%    corr_mat(2,ii)=corr(simulated_series_filtered(y_pos,:,ii)',simulated_series_filtered(c_pos,:,ii)');
%    corr_mat(3,ii)=corr(simulated_series_filtered(y_pos,:,ii)',simulated_series_filtered(i_pos,:,ii)');
%    corr_mat(4,ii)=corr(simulated_series_filtered(y_pos,:,ii)',simulated_series_filtered(k_pos,:,ii)');
%    corr_mat(5,ii)=corr(simulated_series_filtered(y_pos,:,ii)',simulated_series_filtered(h_pos,:,ii)');
%    corr_mat(6,ii)=corr(simulated_series_filtered(y_pos,:,ii)',simulated_series_filtered(productivity_pos,:,ii)');
%end

%Print table with results
%fprintf('\n%-40s \n',title_string)
%fprintf('%-20s \t %11s \t %11s \n','','std(x)','corr(y,x)')
%for ii=1:size(corr_mat,1)
%    fprintf('%-20s \t %3.2f (%3.2f) \t %3.2f (%3.2f) \n',var_names{ii,:},mean(std_mat(ii,:,:),3),std(std_mat(ii,:,:),0,3),mean(corr_mat(ii,:),2),std(corr_mat(ii,:),0,2))
%end
